#############################################################################
# Non linear temperature measures: gdd and kdd
#
#############################################################################

rm(list=ls())
graphics.off()

library(data.table)
library(ggplot2)
library(haven)
library(fixest)
library(xtable)
library(Hmisc)
library("RColorBrewer")

memory.limit(size=500000)

save_dta = "Data_new"
time0 = Sys.time()

#############################################################################
# 1. Data
#############################################################################

load(paste0(save_dta,"/Data_baseline.RData"))

# Bring data with all measures
dt_all = data.table(read_dta("Data_ori/Schlenker/Weather_data_1950_2019.dta"))

# Merge
dim(dt)
dt= merge(dt,dt_all[,.SD,.SDcols=c("fips","year","dday29C","dday0C")],by=c("fips","year"),all.x=TRUE)
dim(dt)
summary(dt[,.(ddr,dday0C,dday29C)])
dim(dt)

# Generate gdd (0-29) and kdd (29+) variables
dt[,gdd:=(dday0C-dday29C)/183]
dt[,kdd:=dday29C/183]
summary(dt[,.(ddr,gdd,kdd)])

save_res = "Results/"
source("Code/0. CODE Auxiliary.R")

#############################################################################
# 2. FE with heterogeneous beta, no gamma (baseline)
#############################################################################

# Estimate
f2a = feols(yield ~  precr + ddr:factor(fips) | factor(fips) + factor(year)^factor(state),data=dt)
summary(f2a)

# Check number of coefficients: 2253 slopes + 1 precr 
nobs(f2a)
dim(f2a$coeftable)
summary(f2a$coeftable)

# Save coefficients
dtc_2a = slopes_dt(f=f2a,names_var="")
dtc_2a[,gamma:=0]
summary(dtc_2a)

#############################################################################
# 3. Specification with heterogeneous gdd and kdd
#############################################################################

# Estimate
f3d = feols(yield ~  precr + precr_sq + gdd:factor(fips) + kdd:factor(fips) | factor(fips) + factor(year)^factor(state) ,data=dt,cluster="state")
summary(f3d)

# Check number of coefficients: 2253*2 beta and gamma + 2 (precr, precr_sq)
nobs(f3d)
2253*2+2
dim(f3d$coeftable)

# Save gdd coefficients in table
nn = rownames(f3d$coeftable)
slopes_pos = grep(":kdd",nn)
dtc_3d = data.table(f3d$coeftable[slopes_pos,])
nn_pos = substr(nn[slopes_pos],13,17)
nn_pos = sub(pattern=":",replacement="",x=nn_pos)
nn_pos = as.numeric(nn_pos)
length(unique(nn_pos))
dtc_3d[,fips:=nn_pos]
setnames(dtc_3d,"Estimate","gamma")
dtc_3d= dtc_3d[,.(gamma,fips)]
summary(dtc_3d)

# Save gdd coefficients in table
nn = rownames(f3d$coeftable)
slopes_pos = grep("gdd:",nn)
dtc_aux = data.table(f3d$coeftable[slopes_pos,])
nn_pos = substr(nn[slopes_pos],17,21)
nn_pos = sub(pattern=":",replacement="",x=nn_pos)
nn_pos = as.numeric(nn_pos)
length(unique(nn_pos))
dtc_aux[,fips:=nn_pos]
setnames(dtc_aux,"Estimate","beta")
dtc_aux= dtc_aux[,.(beta,fips)]
dtc_3d = merge(dtc_3d,dtc_aux,by="fips",all=TRUE)
summary(dtc_3d)
dim(dtc_3d)

rm(f2a,f3d)
time1=Sys.time()
save.image("Data_new/FE_nonlinear.RData")

#############################################################################
# 4. Counterfactual effects of increasing 1C of temperature per day
#############################################################################

rm(list=ls())
load("Data_new/FE_nonlinear.RData")

### Effects for the model with dd (without distinguishing between gdd and kdd)

dt_dd = merge(dt,dtc_2a[,.(fips,coef)],by="fips",all=TRUE)
dim(dt_dd)
setnames(dt_dd,"coef","mg_base")
summary(dt_dd)

#### Effects for model with gdd and kdd

# Bring data of measures in the conterfactuals
dt_1c = data.table(read_dta("Data_ori/Schlenker/Weather_data_1950_2019_1C.dta"))
dt_1c[,gdd1:=(dday0C-dday29C)/183]
dt_1c[,kdd1:=dday29C/183]
summary(dt_1c[,.(gdd1,kdd1)])
dim(dt_1c)
summary(dt_1c)

# Merge two data sets to have original temperature
dt_kdd = dt[,.(fips,year,gdd,kdd)]
summary(dt_kdd)
dt_kdd = merge(dt_kdd,dt_1c[,.(fips,year,gdd1,kdd1)],by=c("fips","year"),all.x=TRUE)
summary(dt_kdd)
summary(dt_kdd[,gdd1-gdd])
summary(dt_kdd[,kdd1-kdd])

# Generate effects
dt_kdd = merge(dt_kdd,dtc_3d,by="fips",all=TRUE)
dt_kdd[,mg_nl:=beta*(gdd1-gdd)+gamma*(kdd1-kdd)]
summary(dt_kdd)
dt_kdd[,c("beta","gamma"):=NULL]
dt_kdd[,c("gdd","kdd"):=NULL]

##### Put together

dim(dt_dd)
dim(dt_kdd)
dt_1c = merge(dt_dd,dt_kdd,by=c("fips","year"),all=TRUE)
dim(dt_1c)
summary(dt_1c)

########################
# Tabulate temperatures
########################

dt_1c[,ddr_1c:=ddr +1]
vnames = c("ddr","gdd","kdd","ddr_1c","gdd1","kdd1")
s=as.matrix(dt_1c[,lapply(.SD,wtd.quantile,probs=c(0.1,0.25,0.5,0.75,0.9),weights=corn_area),
                  .SDcols=vnames],weights=corn_area)
s=rbind(s,as.matrix(dt_1c[,lapply(.SD,wtd.mean,weights=corn_area),.SDcols=vnames]))
s=rbind(s,as.matrix(dt_1c[,lapply(.SD,wtd.var,weights=corn_area),.SDcols=vnames]))
rownames(s)=c(paste0("Percentile ",c(10,25,50,75,90)),"Mean","Variance")
s
l=length(vnames)
dig=rep(3,l+1)
print(xtable(s,type="latex",digits=dig,display=c("s",rep("f",l))),hline.after = NULL,
      file=paste0(save_res,"Table_temp_1c.tex"),include.rownames = TRUE,include.colnames=FALSE,
      sanitize.text.function = function(x){x},only.contents = TRUE)


########################
# Tabulate effects
########################

vnames = c("mg_base","mg_nl")
s=as.matrix(dt_1c[,lapply(.SD,wtd.quantile,probs=c(0.1,0.25,0.5,0.75,0.9),weights=corn_area),
                  .SDcols=vnames],weights=corn_area)
s=rbind(s,as.matrix(dt_1c[,lapply(.SD,wtd.mean,weights=corn_area),.SDcols=vnames]))
s=rbind(s,as.matrix(dt_1c[,lapply(.SD,wtd.var,weights=corn_area),.SDcols=vnames]))
rownames(s)=c(paste0("Percentile ",c(10,25,50,75,90)),"Mean","Variance")
s
l=length(vnames)
dig=rep(3,l+1)
print(xtable(s,type="latex",digits=dig,display=c("s",rep("f",l))),hline.after = NULL,
      file=paste0(save_res,"Table_nonlinear_1c_sub.tex"),include.rownames = TRUE,include.colnames=FALSE,
      sanitize.text.function = function(x){x},only.contents = TRUE)

########################
# Plotting Densities
########################
# Reshape for plotting
summary(dt_1c[,.SD,.SDcols=vnames])
dt_plot = reshape(dt_1c[,.SD,.SDcols=c("fips","year","corn_area",vnames)],direction="long",idvar=c("fips","year"),timevar="G",varying=vnames,sep="_")
table(dt_plot$G)

# Plotting densities: only baseline and fully heterogeneous GDD and KDD
yl=c(0,0.16)
xl=c(-30,10)
gg = ggplot(data=dt_plot)  + theme_bw() +
  theme(title=element_text(size=14),plot.title=element_text(hjust=0.5),
        plot.subtitle=element_text(hjust=0.5), 
        axis.text = element_text(size=18),axis.title = element_text(size=18),
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        legend.position = 'none',legend.text = element_text(size=10)) +
  xlab("Temperature mg effect") + ylab("Density") + xlim(xl) + ylim(yl) +
  scale_color_manual(name = "", breaks = vnames,labels=vnames,
                     values = c("base"="black","ml"="black")) +
  scale_linetype_manual(name = "", breaks = vnames,labels=vnames,
                        values = c("base"="dashed","nl"="solid"))
gg + geom_density(aes(x=mg,weight=corn_area,color=factor(G),linetype=factor(G)),adjust=1.5) 
ggsave(paste0(save_res,"Density_1c_sub.pdf"),width=6.3,height=3.9,units="in")

########################
# Maps 
########################

# Bring US map
load("Data_ori/Data_us_map.RData")
source("Code/0. CODE Auxiliary.R")

# Average of effects per county of increasing 1C
dt_avg = dt_1c[,lapply(.SD,wtd.mean,weights=corn_area),by=fips,.SDcols=vnames]
dim(dt_avg)
summary(dt_avg)

# Effects
qq= quantile(dt_avg$mg_nl,probs=seq(0.05,0.95,0.05))
qq
qq= quantile(dt_avg$mg_base,probs=seq(0.05,0.95,0.05))
qq
bmin= qq[2]
bmax= qq[18]
plot_map(data=dt_avg[,.(fips,mg_base)],trunc=1,bmin=bmin,bmax=bmax,save_name="Map_1c_baseline",palette="OrRd",rev=1,nc=9)
plot_map(data=dt_avg[,.(fips,mg_nl)],trunc=1,bmin=bmin,bmax=bmax,save_name="Map_1c_gdd_het_kdd_het",palette="OrRd",rev=1,nc=9)

# KDD coefficients
dtc_3d[,gammad:=gamma/10]
summary(dtc_3d)
qq = quantile(dtc_3d$gammad,probs = seq(0.05,0.95,0.05))
qq
bmin= qq[2]
bmax= qq[18]
plot_map(data=dtc_3d[,.(fips,gammad)],trunc=1,bmin=bmin,bmax=bmax,save_name="Map_gamma_gdd_het_kdd_het",llabel="Coefficient (in 10s)    ",palette="OrRd",rev=1,nc=9)

# GDD coefficients
summary(dtc_3d)
qq = quantile(dtc_3d$beta,probs = seq(0.05,0.95,0.05))
qq
bmin= qq[2]
bmax= qq[18]
plot_map(data=dtc_3d[,.(fips,beta)],trunc=1,bmin=bmin,bmax=bmax,save_name="Map_beta_gdd_het_kdd_het",llabel="Coefficient    ",palette="OrRd",rev=1,nc=9)


